\documentclass[UTF8]{ctexart}

\usepackage{tcolorbox}

\newtcolorbox{notebox}{
  colback=yellow!10!white,
  colframe=yellow!50!black,
  fonttitle=\bfseries,
  title=注意}

\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhf{}
\rhead{\rightmark}
\cfoot{\thepage}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb} 
\usepackage{amsthm} 
\usepackage{tikz}
\usepackage{xeCJK}
\usepackage{ctex}
\usepackage{url}

\begin{document}
\title{项目作业}
\author{黄紫萱\quad 3220103649}
\date{2023年7月9日}
\maketitle
\newpage

\section{介绍}

\textbf{\large{有限元方法的基本设置}}\\
这是我们实际使用有限元来计算某些东西的第一个示例。我们将求解一个边界值为零但右侧为非零的泊松方程的简单版本：\cite{wangzhan}\\
\begin{align*} -\Delta u &= f \qquad\qquad & \text{在}\ \Omega, \\ u &= 0 \qquad\qquad & \text{在}\ \partial\Omega. \end{align*}
我们将在正方形$\Omega=[−1,1]^2$上解决这个方程，您已经学会了如何在步骤1和步骤2中生成网格。在本程序中，我们只考虑特殊情况$f(x)=1$，并将在下一个教程程序（步骤4）中讨论如何实现更一般的情况。\\
如果您已经了解了有限元方法的基础知识，您将记得我们需要采取的步骤来近似解u通过有限维近似。具体来说，我们首先需要推导出上述方程的弱形式，我们通过将方程乘以测试函数来获得$\varphi$从左边（我们将回到从左边乘法而不是从右边乘法的原因）并在域上积分$\Omega$:
\begin{equation}
-\int_{\Omega} \varphi \Delta u = \int_{\Omega} \varphi f
\end{equation}
这可以通过部件进行集成：
\begin{equation}
\int_{\Omega} \nabla \varphi \cdot \nabla u - \int_{\partial \Omega} \varphi n \cdot \nabla u = \int_{\Omega} \varphi f
\end{equation}
测试功能$\varphi$必须满足相同类型的边界条件（用数学术语来说：它需要来自我们寻求解的集合的切空间），所以在边界上$\varphi=0$因此，我们正在寻找的弱形式读取
\begin{center}
$(\nabla \varphi, \nabla u) = (\varphi, f)$
\end{center}
我们使用了通用符号的地方$(a,b) = \int_{\Omega} a b$然后问题要求一个函数u此语句适用于所有测试函数$\varphi$从适当的空间（这里是空间$H^1$).\\
当然，在一般情况下，我们无法在计算机上找到这样的功能，而是寻求近似值$u_h(x) = \sum_{j} U_j \varphi_j(x)$其中$U_j$是我们需要确定的未知膨胀系数（这个问题的“自由度”），以及$\varphi_i(x)$是我们将使用的有限元形状函数。要定义这些形状函数，我们需要以下内容：\\
\begin{itemize}
  \item 用于定义形状函数的网格。您已经了解了如何在步骤 1 和步骤 2 中生成和操作描述网格的对象。\\
  \item 一个有限元，用于描述我们要在参考单元格上使用的形状函数。\\
  \item 一个 DoFHandler 对象，它枚举网格上的所有自由度，以有限元对象提供的参考单元格描述为基础。您还已经在步骤 2 中了解了如何执行此操作。\\
  \item 一种映射，告知如何从参考单元上的有限元类定义的形状函数中获取实际单元上的形状函数。默认情况下，除非您另有明确说明，否则 deal.II 将为此使用（双向、三）线性映射，因此在大多数情况下，您不必担心此步骤。\\
\end{itemize}

通过这些步骤，我们现在有了一组函数 $\varphi_i$，我们可以定义离散问题的弱形式：寻找一个函数 $u_h$，即寻找上述提到的扩展系数 $U_j$，使得：\\
\begin{center}
$(\nabla \varphi_i, \nabla u_h) = (\varphi_i, f), \quad i = 0, \ldots, N-1.$
\end{center}
请注意，我们在这里遵循从零开始计数的约定，这在 C 和 C++ 中很常见。如果您插入表示 uh(x) 的表示形式 $uh(x) = \sum_{j} U_j \varphi_j(x)$，然后观察到\\
\begin{center}
$(\nabla \varphi_i, \nabla u_h) = (\nabla \varphi_i, \nabla \left[\sum_{j} U_j \varphi_j\right]) = \sum_{j} (\nabla \varphi_i, \nabla [U_j \varphi_j]) = \sum_{j} (\nabla \varphi_i, \nabla \varphi_j) U_j.$
\end{center}
有了这个，问题如下：找到一个向量U因此
$AU=F$
矩阵 A 和右手边项 F 的定义如下：
$A_{ij}= (\nabla \varphi_i, \nabla \varphi_j)$\\
$F_i = (\varphi_i, f)$\\

\textbf{\large{我们应该从左乘还是从右乘以测试函数？}}\\
在我们继续描述如何计算这些量之前，请注意，如果我们从右边乘以测试函数而不是从左边乘以测试函数，那么我们将得到一个形式的线性系统：
\[
U^T A = F^T
\]
使用行向量$F^T$。通过转置这个系统，这当然相当于解决
\[
A^T U = F
\]
这里的情况与上述情况相同，因为 $A=A^T$。但一般情况下，并非如此。为了避免任何混淆，经验表明，习惯上将方程从左边而不是从右边进行乘法（这在数学文献中经常这样做），可以避免一类常见错误，因为矩阵在进行理论和实现比较时不需要进行转置。在本教程的第9步中，我们有一个非对称的双线性形式，从右边还是从左边进行乘法会产生不同的结果。请参考该教程的第9步作为一个示例。\\

\textbf{\large{组装矩阵和右侧向量}}\\
现在我们知道我们需要的是什么（即：持有矩阵和向量的对象，以及计算 $A_{ij}$ 和 $F_i$ 的方法），我们可以看看如何实现它：
\begin{itemize}

\item A 的对象类型是 SparseMatrix，而 U 和 F 的对象类型是 Vector。我们将在下面的程序中看到用于求解线性系统的类。\\
\item 我们需要一种形成积分的方法。在有限元方法中，这通常是使用数值积分来完成的，即积分被一组积分点上的加权和所代替。也就是说，我们首先将关于 Ω 的积分分割成所有单元上的积分，\\
\begin{align*} 
A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j) = \sum_{K \in {\mathbb T}} \int_K \nabla\varphi_i \cdot \nabla \varphi_j, \\ 
F_i &= (\varphi_i, f) = \sum_{K \in {\mathbb T}} \int_K \varphi_i f, 
\end{align*}
然后用正交近似每个单元格的贡献：\\
\begin{align*} 
A^K_{ij} &= \int_K \nabla\varphi_i \cdot \nabla \varphi_j \approx \sum_q \nabla\varphi_i(\mathbf x^K_q) \cdot \nabla \varphi_j(\mathbf x^K_q) w_q^K, \\ 
F^K_i &= \int_K \varphi_i f \approx \sum_q \varphi_i(\mathbf x^K_q) f(\mathbf x^K_q) w^K_q, 
\end{align*}
其中$T \approx \Omega$是近似域的三角测量，$x_q^K$是单元格 K 上的第 q 个正交点，$w_q^K$是第 q 个正交权重。这样做需要什么有不同的部分，我们接下来将依次讨论它们。\\
\item 首先，我们需要一种描述积分点的位置 $x_{Kq}$和权重 $w_{Kq}$ 的方法。它们通常是从参考单元映射而来，类似于形函数的方式，即隐式地使用 MappingQ1 类，或者如果您明确指定，可以通过从 Mapping 类派生的其他类进行映射。在参考单元上描述位置和权重的对象是从 Quadrature 基类派生的。通常，我们会选择一个积分公式（即一组积分点和权重），以便积分与矩阵中的积分完全相等；这可以通过高斯积分公式实现，该公式在 QGauss 类中实现。\\
\item 然后，我们需要一种方法来在单元 K 上评估 $\varphi_i(x_{Kq})$。这就是 FEValues 类所做的工作：它接受一个有限元对象来描述参考单元上的 $\varphi$，一个积分对象来描述积分点和权重，以及一个映射对象（或隐式地接受\textbf{MappingQ1}类），并在位于单元 K 上的积分点上提供形函数及其导数的值，以及进行积分所需的所有其他信息。\\
\end{itemize}
\textbf{FEValues}真是汇集过程中的核心类。你可以将它视为以下方面：\textbf{FiniteElement}和派生类描述形状函数，也就是，无限维对象：函数在每个点都有值。我们需要这个，因为我们希望用函数上的积分进行分析。然而，对于计算机来说，这是一个非常困难的概念，因为它们通常只能处理有限量的信息，所以我们用映射（映射对象）通过在参考单元上定义的点（求积对象）替换积分到实际单元上的点。本质上，我们将问题简化为只需要有限量的信息，即形状函数的值和导数，求积权重，法向量等，仅在有限集合的点上。\textbf{FEValues}类是将这三个组件结合在一起并在特定单元 $K$ 上提供这个有限集合的信息的类。当我们在下面组装线性系统时，你会看到它的作用。

值得注意的是，如果你在应用程序中自己创建了这三个对象，并自己处理信息，也可以达到所有这些目标。然而，这既不更简单（{\textbf{FEValues}类提供了你实际需要的信息），也不更快：\textbf{FEValues}类是高度优化的，只在每个单元上计算你需要的特定信息；如果有任何东西可以从前一个单元重用，那它就会这样做，而且在这个类中有很多代码确保了只要有利，就会缓存这些东西。

这个介绍的最后一部分是提及，获得线性系统后，使用迭代求解器解决，然后进行后处理：我们使用{\textbf{DataOut}类创建一个输出文件，然后可以使用常见的可视化程序进行可视化。
\begin{notebox}
注意：上述所有有限元实现的重要步骤在deal.II中都有对应：该库自然可以分为一些“模块”，覆盖了刚才概述的基本概念。你可以通过本页顶部的标签页访问这些模块。在deal.II手册的首页上也可以看到最基本的概念群的概述。
\end{notebox}
\textbf{\large{关于实现}}\\
尽管这是使用有限元方法可以求解的最简单的方程，但该程序显示了大多数有限元程序的基本结构，并且还可以作为几乎所有后续程序基本上遵循的模板。具体来说，该程序的主类如下所示：
\begin{verbatim}
class Step3
{
  public:
    Step3 ();
    void run ();
 
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
 
    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;
 
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{verbatim}
这遵循了数据封装的面向对象编程口号，即我们尽最大努力将此类的几乎所有内部细节隐藏在外部无法访问的私有成员中。
让我们从成员变量开始：它们遵循我们在上面项目符号中概述的构建块，即我们需要一个三角测量和一个 DoFHandler 对象，以及一个描述我们要使用的形状函数类型的有限元对象。第二组对象与线性代数有关：系统矩阵和右侧以及解向量，以及描述矩阵稀疏模式的对象。这就是此类的全部需求（也是任何稳态偏微分方程求解器所需的基本要素），并且需要在整个程序中生存。与此相反，我们组装所需的 FEValues 对象在整个程序集中都是必需的，因此我们在执行此操作的函数中将其创建为本地对象，并在其末尾再次销毁它。
其次，让我们看一下成员函数。这些也已经形成了几乎所有后续教程程序都将使用的通用结构：

\begin{itemize}
  \item \textbf{make\_grid()}：这就是可以称为预处理函数的函数。顾名思义，它设置存储三角测量的对象。在后面的示例中，它还可以处理边界条件、几何形状等。
  \item \textbf{setup\_system()}：这是设置解决问题所需的所有其他数据结构的功能。特别是，它将初始化 DoFHandler 对象并正确调整与线性代数相关的各种对象的大小。此函数通常与上面的预处理函数分开，因为在与时间相关的程序中，只要网格被自适应细化，就可以至少每隔几个时间步调用它（我们将在步骤 6 中看到如何操作）。另一方面，在上面的预处理函数中设置网格本身只在程序开始时完成一次，因此被分离到它自己的函数中。
  \item \textbf{assemble\_system()}：然后，这就是计算矩阵和右手边的内容的地方，如上面的介绍所述。因为用这个线性系统做事在概念上与计算它的条目非常不同，所以我们将其从下面的函数中分离出来。
  \item \textbf{solve()}：然后，这就是我们计算线性系统$AU=F$的解$U$的函数。在当前的程序中，这是一个简单的任务，因为矩阵很简单，但是当问题不再那么平凡时（例如，在你对库有了更多的了解后，看看step-20、step-22或step-31），它将成为程序大小的重要部分。
  \item \textbf{output\_results()}：最后，当你计算出一个解决方案后，你可能想要做点什么。例如，你可能想要以可以可视化的格式输出它，或者你可能想要计算你感兴趣的数量：比如，换热器中的热通量，翼的空气摩擦系数，最大的桥梁负载，或者仅仅是一个点的数值解的值。因此，这个函数就是用于后处理你的解的地方。
\end{itemize}
所有这些都由单一的公共函数（除了构造函数）run()函数维系在一起。它是从创建这种类型的对象的地方被调用的函数，它是按照正确的顺序调用所有其他函数的函数。将这个操作封装到run()函数中，而不是从main()中调用所有其他函数，确保你可以改变这个类内部职责分离的实现方式。例如，如果其中一个函数变得太大，你可以将它分成两个，你必须关心的改变的唯一的地方就在这个类内部，而不是在其他任何地方。

如上所述，您将在以下许多教程程序中再次看到这种一般结构 - 有时在函数名称的拼写方面有变体，但本质上是功能分离的顺序。\\

\textbf{\large{关于类型的说明}}\\

deal.II通过命名空间类型中的别名在{\textbf{types}命名空间中定义了许多整型类型。在上一句中，“积分”一词用作与名词“整数”对应的形容词。它不应与表示曲线或表面下的面积或体积的名词“积分”混淆。形容词“积分”在C++世界中广泛用于“积分类型”、“积分常数”等上下文中.特别地，在这个程序中，你会在几个地方看到\textbf{types::global\_dof\_index}：一个用于表示自由度全局索引的整数类型，即，定义在三角剖分顶部的DoFHandler对象中特定自由度的索引（与特定单元内特定自由度的索引相对）。对于当前程序（以及几乎所有的教程程序），你全局上会有几千到可能几百万个未知数（对于$Q_1$元素，你在2d的每个单元上会有4个，3d中会有8个）。因此，一个允许存储足够大的全局DoF索引数的数据类型是无符号整数，因为它允许存储0到略大于40亿的数字（在大多数系统中，整数是32位的）。实际上，这就是\\textbf{types::global\_dof\_index}。

那么，为什么不立即使用呢？deal.II在7.3版本之前是这样做的。然而，deal.II支持非常大的计算（通过在step-40中讨论的框架），当跨越几千个处理器时，可能会有超过40亿个未知数。因此，有些情况下，无符号整数不够大，我们需要一个64位的无符号整数类型。为了实现这一点，我们引入了\textbf{types::global\_dof\_index}，它默认定义为简单的无符号整数，而如果需要，可以通过在配置期间传递特定的标志将其定义为无符号长整数（参见ReadMe文件）。

这涵盖了技术方面。但也有文档的目的：在库和构建在其上的代码中的任何地方，如果你看到一个地方使用数据类型{types::global\_dof\_index}，你立即知道被引用的量实际上是一个全局dof索引。如果我们只是使用无符号整数（也可能是一个本地索引，一个边界指示器，一个材料id等），则不会有这样的含义。立即知道一个变量指的是什么也有助于避免错误：如果你看到一个\textbf{types::global\_dof\_index}类型的对象被赋值给\textbf{types::subdomain\_}\\
\textbf{id}类型的变量，即使它们都由无符号整数表示，编译器也不会抱怨，那么很明显，一定有一个错误。

在更实际的意义上,该类型的存在意味着在组装过程中,我们创建一个$4\times 4$矩阵(在2d中,使用$Q_1$元件),表示我们目前所在单元的贡献,然后我们需要将这个矩阵的元素加到全局(系统)矩阵的适当元素上。为此,我们需要获得与当前单元本地相关的自由度的全局索引,对此我们将始终使用以下代码片段:

\begin{verbatim} 
cell->get_dof_indices (local_dof_indices);  
\end{verbatim}

其中local\_dof\_indices声明为

\begin{verbatim}
std::vector<types::global_dof_index> local_dof_indices (fe.dofs_per_cell);
\end{verbatim}

这个变量的名称可能有点误导——它代表“那些本地在当前单元上定义的自由度的全局索引”——但是整个库中统一地将保存此类信息的变量命名为这样。

\begin{notebox}
\textbf{types::global\_dof\_index} 不仅是这个命名空间中定义的唯一类型。相反,这里有一整套类型,包括 
\textbf{types::subdomain\_id}、\textbf{types::boundary\_id} 和 \textbf{types::material\_id}。所有这些都是整数数据类型的别名,但如上所述,它们在整个库中使用,以便(i)变量的意图更容易辨识,和(ii)如果需要,可以将实际类型更改为更大的类型,而无需遍历整个库并确定unsigned int的特定用途是否对应于材料指示器等。
\end{notebox}

\section*{问题和代码的描述}
\indent 此程序旨在解决二维拉普拉斯方程。它涉及以下主要步骤：网格生成、系统设置、系统组装、方程求解和结果输出。\\
\indent 代码也主要分为5个环节：\\
\indent 网格生成（make grid）：首先创建了一个位于 [-1,1]范围内的立方体网格，并进行了五次全局细化。全局细化是将每个单元格细分为更小单元格的过程。最后，打印活动单元格的数量。\\
\indent 系统设置（setup system）：首先为网格上的自由度（DoF）分配空间。然后生成稀疏模式以存储即将生成的稀疏矩阵的结构。最后，初始化系统矩阵、解向量和右侧向量。\\
\indent 系统组装(assemble system)：计算系统矩阵和右侧向量。系统矩阵是通过循环遍历每个单元格并计算单元格矩阵和右侧向量，然后将它们添加到全局矩阵和全局右侧向量中得到的。同时，应用边界条件。\\
\indent 方程求解（solve）：使用共轭梯度法和简单的恒等预处理器来解决线性系统。\\
\indent 结果输出(output results)：最后，这部分将解输出到一个gnuplot文件中，以便于后续的可视化。\\
\indent main函数部分先设置了日志的深度，然后创建了一个Step3类的实例并运行run方法，开始执行上述步骤。\\

\textbf{\Large{生成图片如下：}}\\
\includegraphics{hzx.eps}


\bibliographystyle{plain}
\bibliography{wenxian.bib}
\end{document}
